function solvestep(CFL,n,opt,method)
% SOLVESTEP - 用TVD格式计算台阶问题
%   CFL:为CFL系数。
%   n: 为纵向网格个数。h=1/n;为了方便边界的处理要求n必须为5的倍数。
%   opt:debug选项 
%      1:debug 模式，画出计算过程的中的rho。
%      else: 只画出最后的rho。
%   method:方法选项
%      1：ULT1(default)
%      2：ULT1C
    if nargin <4
        method =1;
    end
    
    if opt == 1
        pause on
    end
    h=1/n;
    m = 3*n;
    global U;
    global t;
    xr = m/5;                           % 台阶区域为[xr,m]*[0,yl];
    yl = n/5;                           % 
    U= cell(n,m);
    
    % 初始化。
    u = Variables2U(1.4,3,0,1);
    %u =[1.4;3;0;6.55];
    for i =1:n
        for j = 1:m
            if i>yl || j<=xr
                U{i,j}=u;
            else
                U{i,j}=nan*ones(4,1);
            end
        end
    end
    t =0;
    T_final = 4;
    dots = 0:T_final/50:T_final;
    dot = 1;
    fprintf('%5.1f%%',0);
    tau2 = nan;
    while(t < T_final)
        if(t>dots(dot))
            PrintPercentage(t/T_final);
            dot=dot+1;
        end
        eigen = maxeigen();
        temp = min(T_final-t,h/eigen);
        if isnan(tau2)
            tau1=CFL*temp/2;
        else
            tau1 = CFL*temp-tau2;
        end
        lambda = tau1/h;
        t = t+tau1;        
        %Lx x方向做一次TVD迭代。
        XTVD(lambda,method);
        modify
        togglestep
        
        %Ly y方向做一次TVD迭代
        YTVD(lambda,method);
        modify
        togglestep
        if opt == 1
            showrho
        end
        if t>= 4
           break
        end
        eigen = maxeigen();
        temp = min(T_final-t,h/eigen);
        tau2 = CFL*temp-tau1;
        lambda = tau2/h;
        t = t+tau2;
        %Ly again
        YTVD(lambda,method);
        modify
        togglestep
        
        %Lx again
        XTVD(lambda,method);
        modify
        togglestep
        if opt == 1
            showrho
        end
    end
    PrintPercentage(1);
    fprintf('\n');
    showrho()
    h = gcf;
    hgexport(h,['n=',num2str(n)]);
end
function PrintPercentage(per)
% PRINTPERCENTAGE - 
%   
    fprintf('\b\b\b\b\b\b');
    fprintf('.');
    per = 100*per;
    fprintf('%5.1f%%',per);    
end

function XTVD(lambda,method)
% XTVD - 
%   
    global U;
    [n,m] = size(U);
    % 台阶下
    for i =1:n/5                         
        U(i,1:m/5) = LTVD(U(i,1:m/5),2,lambda,method);
    end
    % 台阶上面
    for i = (n/5+1):n                   
        U(i,:) = LTVD(U(i,:),1,lambda,method);
    end
end
function YTVD(lambda,method)
% YTVD - 
%   
    global U;
    [n,m] = size(U);
    % 台阶左边。
    for i = 1:m/5                        
        Ut= LTVD(U(:,i),3,lambda,method);
        U(:,i)=Ut(:);
    end
    % 台阶右边。
    for i = (m/5+1):m                   
        Ut = LTVD(U((n/5+1):n,i),3,lambda,method);
        U((n/5+1):n,i)=Ut(:);
    end 
    
end

function showrho()
% SHOWRHO - 
%   
    global U;
    global t;
    [n,m]=size(U);
    h = 1/n;
    rho=zeros(n,m);
    u = zeros(n,m);
    v=zeros(n,m);
    p=zeros(n,m);
    for i = 1:n
        for j = 1:m
            [rho(i,j),u(i,j),v(i,j),p(i,j)]=U2Variables(U{i,j});
        end
    end
    x = h/2:h:3-h/2;
    y = h/2:h:1-h/2;
    if n<16 
        %myprint(rho)
    end
    subplot(2,1, 1)
    contourf(x,y,rho,20);
    hold on; set(gca,'DataAspectRatio',[1 1 1]);set(gca,'Xtick',0:0.2:3);
    title(['Density t=',num2str(t)]);
    xlabel('x');ylabel('y');

    %h = gcf;
    %hgexport(h,'Density');
    hold off
    subplot(2, 1, 2)
    contour(x,y,rho,20);
    hold on; set(gca,'DataAspectRatio',[1 1 1]);set(gca,'Xtick',0:0.2:3);
    title(['Density t=',num2str(t)]);
    xlabel('x');ylabel('y');

    hold off
    %pause
    pause(0.01)
end
function myprint(M)
% MYPRINT - 
%   
    [n,m]=size(M);
    for i=1:n
        for j=1:m
            fprintf('%4.4e ',M(i,j) );
        end
        fprintf('\n');
    end    
    fprintf('\n')
end

function togglestep()
% TOGGLESTEP - 处理台阶角。
%   
    global U;
    [n,m] = size(U);
    %处理台阶角：
    yl = n/5+1;
    xr = m/5+1;
    b=[yl,xr;
       yl,xr+1;
       yl,xr+2;
       yl,xr+3;
       yl+1,xr;
       yl+1,xr+1;
      ];
    [rhoa,ua,va,pa]=U2Variables(U{n/5,m/5});
    for i=1:6
        [rhob,ub,vb,pb] = U2Variables(U{b(i,1),b(i,2)});
        if pa~=0
            rhob =rhoa*(pb/pa)^(1/1.4);
        end
        s1=(ua^2+va^2)/2;
        s2=(ub^2+vb^2)/2;
        if s2 ~=0
            alph =(s1+1.4*(pa/rhoa-pb/rhob))/s2;
        else
            alph =1;
        end
        U{b(i,1),b(i,2)}=[ rhob,alph*rhob*ub,alph*rhob*vb,2.5*pb+alph^2*rhob*s2]';
    end
       
end

function modify()
% MODIFY - 修正 负的压强和密度。
%   
    global U;
    [n,m] = size(U);
    sp = 1e-8;
    srho = 1e-8;
    for i = 1:n
        for j=1:m
            if i<=n/5 && j>m/5
                continue
            end
              
            [rho,u,v,p]=U2Variables(U{i,j});
            if sum(isnan([rho,u,v,p]))
                [rho,u,v,p]
                [i,j]
                showrho
                pause
            end            
            if rho <0
                rho
                rho = srho;
            end
            if p <0
                p
                p = sp;
            end
            U{i,j}=Variables2U(rho,u,v,p);
        end
    end

end

function eigen = maxeigen()
% MAXEIGEN - 求出整个区域上的最大特征值，以协助计算时间步长。
%   
    global U;
    [n,m]=size(U);
    eigen =0;
    for i=1:n,
        for j=1:m,
            if i > n/5 || j< m/5
                [rho,u,v,p]=U2Variables(U{i,j});
                a = sqrt(1.4*p/rho);
                eigen = max(abs([u,v,u-a,u+a,v-a,v+a,eigen]));
            end
        end
    end
end

